Pyromania—Life after a fiery death: fire and plant biomass loading affect dissolved organic matter in experimental ponds
CJ Spiegel et al. 2023, Global Change Biology
11/06/2023
Project overview
The goal of the Pyromania project is to test how terrestrial subsides (plant biomass loading) and burning influence DOM composition and degradation. We used a manipulative experiment to assess a range of plant material quantities (0-400g per tank) and fire treatment (burned vs unburned material) and the non-linearity of these effects on aquatic systems through 4 time-point samplings. We used 400L aquatic experimental ponds and ran the experiment for ~90d in 2021-2022.
Figure. Images showing
(a) controlled burning of plant material,
(b) example of low severity burn (c)
packed plant material in mesh bags, (d) experimental
pond set up, (e) plant material added to ponds, and
(f) in situ incubations.
DATASETS
This dataset is among 3 to be generated for the project and focuses
on:
- Dissolved organic carbon (DOC) and
- Fluorescence and UV-vis spectroscopic absorbance
indices
- DOC Degradation from either photo- or microbial decomposition
- Plant biomass decomposition as % mass loss
TIME POINTS
- Time-0 (T0) = before the addition of plant
materials, plankton were stocked at this stage
- Time-1 (T1) = 10 days after the addition of plant
materials
- Time-2 (T2) = 31 days after …
- Time-3 (T3) = 59 days after …
- Time-4 (T4) = 89 days after …
A Quick note on GAMS
Fitting via restricted maximum likelihood (REML) is advised to give stable results. Fitting the smoothing parameter (sp) to user defined levels not advisable (but certainly is done). Number of basis functions important as well, small # of basis functions (k) will limit wiggliness (and make the model more linear than should be). This wiggliness is set by using the k values. Don’t want too high, or too low…
When using the non-linear smoothing, this is the s(xxxx). Can add a new predictor by adding another smooth function with s(x1) + s(x2). When the variable is inside the smooth function, this is to account for the nonlinear shape. Two smoothing together and you have ADDITIVE Nonlinear modeling. Not all variables need to be non-linear, but adding a predictor outside of the s() then this is allowing for linearity.
In practice, continuous variables are rarely linear in GAMs, automatic smoothing will make a linear shape if setting high smoothing function (ie.. s =1000), then a non-linear variable appears linear. However, setting categorical variables as predicots linear is more common.
Factor-smooth interaction also an option, this is when s(x1 by = fac), this sets different smoothing for different variables of “fac”. Usually, the different smoothers are combined with a varying intercept incase the different categories are different in means and slopes, this would be by adding the s(x1 by = fac) + fac, where the +fac allows for a different slope.
EDF - effective degrees of freedom, equate with wiggliness, with edf =1 as a straight line, and higher edfs as more wiggly. GAM smoother significance is described as not being able to draw a horizontal line through the data. Always check concurvity, which is the collinearity with models from 0-1.
Import DOC data
Import the data for DOC and run a loop to clean up all files. This creates a single dataframe (df) by taking the raw data files, aligning metadata, and filtering to make a clean df for GAMs.
- Analyze DOC at each sampling time point.
- Run model selection and produce plots for each individual timepoint
- Graph the timepoints as a 5 panel figure.
## import treatment IDs
IDs<-read.csv("data/treatment.IDs.csv")
##### grab files in a list
Total.DOC.files <- list.files(path="data/DOC.TN", pattern = "csv$", full.names = T)
Total.DOC.files
## [1] "data/DOC.TN/DOC_T0.csv" "data/DOC.TN/DOC_T1.csv" "data/DOC.TN/DOC_T2.csv"
## [4] "data/DOC.TN/DOC_T3.csv" "data/DOC.TN/DOC_T4.csv"
##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="data/DOC.TN", pattern = "csv$", full.names = F))
file.names
## [1] "DOC_T0" "DOC_T1" "DOC_T2" "DOC_T3" "DOC_T4"
############ formatting all data in for loop
for(i in 1:length(Total.DOC.files))
{
data<-read.csv(Total.DOC.files[i], sep=",")
data<-data[,c(1,3,4)] # removed columns we don't need
data$File<-Total.DOC.files[i]
colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
data$Tank<- IDs$Tank
data$Tank<-as.numeric(as.character(data$Tank)) # make the column of samples all numeric
data <- data[!is.na(as.numeric(as.character(data$Tank))),] # remove all rows that aren't numeric/tanks
data$Treatment<-IDs$Treatment
data$plant.mass..g<-IDs$plant.mass..g
make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
}
########## this is the end of the loop
ls() #see all dfs you've made, the above will be df matching their file names
## [1] "data" "DOC_T0" "DOC_T1" "DOC_T2"
## [5] "DOC_T3" "DOC_T4" "Fig.formatting" "file.names"
## [9] "i" "IDs" "Total.DOC.files"
#Combine files from loop to single df
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)
DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names
# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27)
#alternatively
# remove the 9 letters ('^.) at start
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File)
# if equals DOC_T0_11052021 then, T0, if not then T1
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
ifelse (DOC.df$File=="DOC_T1.csv", "T1",
ifelse (DOC.df$File=="DOC_T2.csv", "T2",
ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))
#rearrange
DOC.df<- DOC.df %>%
select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L)
DOC.df$Treatment<-as.factor(DOC.df$Treatment)DOC analysis
- First, run GAMs to see what is happening and where differences
lie.
- We will then make the five panel figure
######## T0 model
m1.DOC.T0.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.DOC.T0.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.DOC.T0.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.DOC.AIC<-AIC(m1.DOC.T0.gam, m2.DOC.T0.gam, m3.DOC.T0.gam)
DOC.sum_T0=summary(m3.DOC.T0.gam)
DOC.sum_T0=anova.gam(m3.DOC.T0.gam)
gam.check(m3.DOC.T0.gam, rep=1000)
draw(m3.DOC.T0.gam)
concrvity(m3.DOC.T0.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T0.gam, all.terms = TRUE, page=1)
DOC.diff.T0<-plot_difference(
m1.DOC.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
) +
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T0.mod.plot<-
plot_smooths(
model = m3.DOC.T0.gam,
series = plant.mass..g
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-0") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
######## T1 model
m1.DOC.T1.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.DOC.T1.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.DOC.T1.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.DOC.AIC<-AIC(m1.DOC.T1.gam, m2.DOC.T1.gam, m3.DOC.T1.gam)
summary(m1.DOC.T1.gam)
T1.DOC.anova=anova(m1.DOC.T1.gam)
gam.check(m1.DOC.T1.gam, rep=1000)
draw(m1.DOC.T1.gam)
concrvity(m1.DOC.T1.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T1.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T1<-plot_difference(
m1.DOC.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T1.mod.plot<-
plot_smooths(
model = m1.DOC.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-10") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T2
m1.DOC.T2.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.DOC.T2.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.DOC.T2.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.DOC.AIC<-AIC(m1.DOC.T2.gam, m2.DOC.T2.gam, m3.DOC.T2.gam)
summary(m1.DOC.T2.gam)
T2.DOC.anova=anova.gam(m1.DOC.T2.gam)
gam.check(m1.DOC.T2.gam, rep=1000)
draw(m1.DOC.T2.gam)
concrvity(m1.DOC.T2.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T2.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T2<-plot_difference(
m1.DOC.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T2.mod.plot<-
plot_smooths(
model = m1.DOC.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-31") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T3
m1.DOC.T3.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.DOC.T3.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.DOC.T3.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.DOC.AIC<-AIC(m1.DOC.T3.gam, m2.DOC.T3.gam, m3.DOC.T3.gam)
summary(m1.DOC.T3.gam)
T3.DOC.anova=anova.gam(m1.DOC.T3.gam)
gam.check(m1.DOC.T3.gam, rep=1000)
draw(m1.DOC.T3.gam)
concrvity(m1.DOC.T3.gam)
par(mfrow = c(2, 2))
plot(m1.DOC.T3.gam, all.terms = TRUE, page=1)
# model predictions
DOC.diff.T3<-plot_difference(
m1.DOC.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
###########
#plot for the model output on rawdata
DOC.T3.mod.plot<-
plot_smooths(
model = m1.DOC.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-59") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
########## T4
m1.DOC.T4.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.DOC.T4.gam=gam(DOC..mg.L ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.DOC.T4.gam=gam(DOC..mg.L ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.DOC.AIC<-AIC(m1.DOC.T4.gam, m2.DOC.T4.gam, m3.DOC.T4.gam)
summary(m3.DOC.T4.gam)
T4.DOC.anova=anova(m1.DOC.T4.gam)
gam.check(m3.DOC.T4.gam, rep=1000)
draw(m3.DOC.T4.gam)
concrvity(m3.DOC.T4.gam)
par(mfrow = c(2, 2))
plot(m3.DOC.T4.gam, all.terms = TRUE, page=1)
###########
#plot for the model output on rawdata
DOC.T4.mod.plot<-
plot_smooths(
model = m3.DOC.T4.gam,
series = plant.mass..g
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=DOC..mg.L, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim=c(0, 60)) +
ggtitle("Day-89") +
ylab("DOC (mg/L)") +
xlab("plant material (g)") +
Fig.formatting
# model predictions
DOC.diff.T4<-plot_difference(
m1.DOC.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-10,13))+
Fig.formatting
#labels
mod.DOC.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
#adding columns
mod.DOC.column<- data.frame(mod.DOC.rep)
AIC.DOC<-bind_rows(T1.DOC.AIC, T2.DOC.AIC, T3.DOC.AIC, T4.DOC.AIC)
AIC.DOC.mod<-as.data.frame(cbind(mod.DOC.column, AIC.DOC))
names(AIC.DOC.mod)[1]="Model"
AIC.DOC.mod <- AIC.DOC.mod %>%
# Creating an empty column:
add_column(Variable= c("DOC concentration"), .before="Model")
write.csv(AIC.DOC.mod, "output/AIC models/AIC.DOC.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/DOC.mod.AIC.pdf", width = 6.8, height = 3.8) # Export PDF
grid.table(AIC.DOC.mod, rows = NULL)
dev.off()DOC output
This generates a table with GAM summaries
####Parametric output
DOC.para.sums=as.data.frame(rbind(T1.DOC.anova$pTerms.table, T2.DOC.anova$pTerms.table, T3.DOC.anova$pTerms.table, T4.DOC.anova$pTerms.table))
names(DOC.para.sums)[1]="df/edf"
#create empty column for ref.df
DOC.para.sums <- DOC.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
DOC.smooth.sums=as.data.frame(rbind(T1.DOC.anova$s.table, T2.DOC.anova$s.table, T3.DOC.anova$s.table, T4.DOC.anova$s.table))
names(DOC.smooth.sums)[1]="df/edf"
#create empty column for ref.df
DOC.smooth.sums <- DOC.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
DOC.mod.allsums=rbind(DOC.para.sums,DOC.smooth.sums)
rownames(DOC.mod.allsums)<-c(1:12)
#create new ID column
DOC.mod.allsums <- DOC.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("DOC concentration"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/DOC.mod.allsums.pdf", width = 9.3, height = 3.8) # Export PDF
grid.table(DOC.mod.allsums, rows= NULL)
dev.off()DOC model plot
###### compile the plots with effect plots}
extract.legend <- get_legend(
# create some space to the left of the legend
DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))
DOC.alltimes<-plot_grid(
DOC.T0.mod.plot+ theme(legend.position = "none"),
DOC.T1.mod.plot+ theme(legend.position = "none"),
DOC.T2.mod.plot+ theme(legend.position = "none"),
DOC.T3.mod.plot+ theme(legend.position = "none"),
DOC.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,8,3), ncol=6)
ggsave("figures/DOC.alltimes.mod.pdf", height=4, width=18)
DOC.alltimes
Figure. Dissolved organic carbon (DOC) concentration across time in treatments receiving burned and unburned plant material at the start of the experiment and four sampling periods after plant material added. Black lines with gray confidence intervals indicate global smoothers across all data points; solid (burned) and dotted (unburned) black lines together represent treatment-level intercepts with global smoothers; colored lines indicate factor-smooths that vary between treatments.
DOC diff plot
From the GAMs, we can plot exactly where the differences lie.
DOC.mod.diffs<-plot_grid(
DOC.diff.T0+ theme(legend.position = "none")+ ggtitle("Day-0"),
DOC.diff.T1+ theme(legend.position = "none")+ ggtitle("Day-10"),
DOC.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
DOC.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
DOC.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8,8), ncol=5)
DOC.mod.diffs
Figure. Differences in dissolved organic carbon (DOC) concentration for fitted smooth functions between treatments (burned vs. unburned) across plant biomass. Difference in trends shown in red with approximate 95% confidence intervals.
ggsave("figures/DOC.mod.diffs.pdf", height=4, width=16)Import fluorescence and UV-vis absorbance data
Import the data for fluorescence and UV-vis absorbance, and run a loop to clean up all files. This creates a single dataframe (df) by taking the raw data files, aligning metadata, and filtering to make a clean df for GAMs.
- We will analyze fluorescence and UV-vis absorbance across the biomass gradient at each sampling time point.
- Run model selection and produce plots for each index.
- Compile all indices into one plot.
Below are brief descriptions of each index:
SUVA: the specific ultraviolet absorbance at 254 nm (SUVA254) is used as a proxy for DOM aromaticity (Weishaar et al. 2003).
Sr: the ratio of the slope (SR) parameters (S275-295:S350-400) as a proxy for molecular weight (MW) (Helms et al. 2008).
HIX: the humification index (HIX = ratio of areas under the emission curve at 435-480 nm and 300-345 nm plus 435-480 nm at an excitation wavelength 245 nm) that described the degree of humification (high HIX = more terrestrial humic-like material) (Zsolnay et al. 1999).
BIX: the freshness index (BIX), representing the degree to which organic matter has undergone processing and degradation, as the ratio of emission intensity at 380 nm to the maximum emission intensity between 420 and 435 nm at excitation wavelength 310 nm (Parlanti 2000).
FI: the fluorescence index (FI = emission intensity at 470 nm divided with that of 520 nm at 370 nm excitation) which indicated source material (microbial-precursor FI value ≈ 1.8 or terrestrial-precursor FI ≈ 1.3) (McKnight et al. 2001).
Import SUVA data
These files are uploaded separately.
###### SUVA Time point 0
suva.t0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T0.csv")
suva.t0.df$Sample.Name=gsub("\\..*", "", suva.t0.df$Sample.Name)
names(suva.t0.df)[1]="Tank"
suva.t0.df=suva.t0.df[-1,]
suva.t0.df$Tank=gsub("T", "", suva.t0.df$Tank)
suva.t0.df$Tank=as.numeric(suva.t0.df$Tank)
suva.t0.df=suva.t0.df[order(suva.t0.df$Tank),]
#add important info
suva.t0.df$DOC..mg.L=DOC_T0$DOC..mg.L
suva.t0.df$Time.point="T0"
suva.t0.df$Treatment=IDs$Treatment
suva.t0.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t0.df$SUVA=(suva.t0.df[suva.t0.df$Tank,]$abs254/suva.t0.df[suva.t0.df$Tank,]$DOC..mg.L)*100
###### SUVA Time point 1
library(dplyr)
suva.t1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T1.csv")
suva.t1.df=suva.t1.df[-1,]
suva.t1.df=suva.t1.df[-1,]
names(suva.t1.df)[1]="Tank"
suva.t1.df$Tank=gsub("TA", "", suva.t1.df$Tank)
suva.t1.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t1.df$Tank)
suva.t1.df$Tank=sub("\\..*", "", suva.t1.df$Tank)
suva.t1.df$Tank=as.numeric(suva.t1.df$Tank)
suva.t1.df=suva.t1.df[order(suva.t1.df$Tank),]
#add important info
suva.t1.df$DOC..mg.L=DOC_T1$DOC..mg.L
suva.t1.df$Time.point="T1"
suva.t1.df$Treatment=IDs$Treatment
suva.t1.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t1.df$SUVA=(suva.t1.df[suva.t1.df$Tank,]$abs254/suva.t1.df[suva.t1.df$Tank,]$DOC..mg.L)*100
###### SUVA Time point 2
suva.t2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T2_d1A.edit.csv")
names(suva.t2.df)[1]="Tank"
suva.t2.df$Tank=gsub("TA", "", suva.t2.df$Tank)
suva.t2.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t2.df$Tank)
suva.t2.df$Tank=sub("\\..*", "", suva.t2.df$Tank)
suva.t2.df$Tank=as.numeric(suva.t2.df$Tank)
suva.t2.df=suva.t2.df[order(suva.t2.df$Tank),]
#add important info
suva.t2.df$DOC..mg.L=DOC_T2$DOC..mg.L
suva.t2.df$Time.point="T2"
suva.t2.df$Treatment=IDs$Treatment
suva.t2.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t2.df$SUVA=(suva.t2.df[suva.t2.df$Tank,]$abs254/suva.t2.df[suva.t2.df$Tank,]$DOC..mg.L)*100
###### SUVA Time point 3
suva.t3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T3.csv")
names(suva.t3.df)[1]="Tank"
suva.t3.df$Tank=gsub("TA", "", suva.t3.df$Tank)
suva.t3.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t3.df$Tank)
suva.t3.df$Tank=sub("\\..*", "", suva.t3.df$Tank)
suva.t3.df$Tank=as.numeric(suva.t3.df$Tank)
suva.t3.df=suva.t3.df[order(suva.t3.df$Tank),]
#add important info
suva.t3.df$DOC..mg.L=DOC_T3$DOC..mg.L
suva.t3.df$Time.point="T3"
suva.t3.df$Treatment=IDs$Treatment
suva.t3.df$plant.mass..g=IDs$plant.mass..g
class(suva.t3.df$DOC..mg.L)
#calculate SUVA
suva.t3.df$SUVA=(suva.t3.df[suva.t3.df$Tank,]$abs254/suva.t3.df[suva.t3.df$Tank,]$DOC..mg.L)*100
###### SUVA Time point 4
suva.t4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/Results_T4.csv")
names(suva.t4.df)[1]="Tank"
suva.t4.df$Tank=gsub("TA", "", suva.t4.df$Tank)
suva.t4.df$Tank=sub("(.{2})(.*)", "\\1.\\2", suva.t4.df$Tank)
suva.t4.df$Tank=sub("\\..*", "", suva.t4.df$Tank)
suva.t4.df$Tank=as.numeric(suva.t4.df$Tank)
suva.t4.df=suva.t4.df[order(suva.t4.df$Tank),]
#add important info
suva.t4.df$DOC..mg.L=DOC_T4$DOC..mg.L
suva.t4.df$Time.point="T4"
suva.t4.df$Treatment=IDs$Treatment
suva.t4.df$plant.mass..g=IDs$plant.mass..g
#calculate SUVA
suva.t4.df$SUVA=(suva.t4.df[suva.t4.df$Tank,]$abs254/suva.t4.df[suva.t4.df$Tank,]$DOC..mg.L)*100
###### Combine all SUVA time points
suva.df=rbind(suva.t0.df, suva.t1.df, suva.t2.df, suva.t3.df, suva.t4.df)Import the rest of the fluorescence and UV-vis absorbance data
# Time 0
index.T0.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T0_indices.csv")
index.T0.df$Time.point="T0"
names(index.T0.df)[6]="Humification.Index"
index.T0.df=index.T0.df[-c(7:11,13,14)]
#Time 1
index.T1.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T1_indices.csv")
index.T1.df$Time.point="T1"
index.T1.df=index.T1.df[-c(7:11,13,14)]
#Time 2
index.T2.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T2_indices.csv")
index.T2.df$Time.point="T2"
index.T2.df=index.T2.df[-c(7:11,13,14)]
#Time 3
index.T3.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T3_indices.csv")
index.T3.df$Time.point="T3"
index.T3.df=index.T3.df[-c(7:11,13,14)]
#Time 4
index.T4.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/EEMs/T4_indices.csv")
index.T4.df$Time.point="T4"
index.T4.df=index.T4.df[-c(7:11,13,14)]
#combine all index time points
index.df=rbind(index.T0.df, index.T1.df, index.T2.df, index.T3.df, index.T4.df)Add indices to original DOC df
Now we will have all of our data in one single df.
#import and do a loop to clean up all files and make stacked data in single df
DOC_T0$SUVA=suva.df[suva.df$Time.point=="T0",]$SUVA
DOC_T1$SUVA=suva.df[suva.df$Time.point=="T1",]$SUVA
DOC_T2$SUVA=suva.df[suva.df$Time.point=="T2",]$SUVA
DOC_T3$SUVA=suva.df[suva.df$Time.point=="T3",]$SUVA
DOC_T4$SUVA=suva.df[suva.df$Time.point=="T4",]$SUVA
#HIX
DOC_T0$hix=index.df[index.df$Time.point=="T0",]$Humification.Index
DOC_T1$hix=index.df[index.df$Time.point=="T1",]$Humification.Index
DOC_T2$hix=index.df[index.df$Time.point=="T2",]$Humification.Index
DOC_T3$hix=index.df[index.df$Time.point=="T3",]$Humification.Index
DOC_T4$hix=index.df[index.df$Time.point=="T4",]$Humification.Index
#Fluorescence
DOC_T0$fl=index.df[index.df$Time.point=="T0",]$Fluroscence.Index
DOC_T1$fl=index.df[index.df$Time.point=="T1",]$Fluroscence.Index
DOC_T2$fl=index.df[index.df$Time.point=="T2",]$Fluroscence.Index
DOC_T3$fl=index.df[index.df$Time.point=="T3",]$Fluroscence.Index
DOC_T4$fl=index.df[index.df$Time.point=="T4",]$Fluroscence.Index
#slope ratio
DOC_T0$sr=suva.df[suva.df$Time.point=="T0",]$Spectral.Slope.Ratio
DOC_T1$sr=suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio
DOC_T2$sr=suva.df[suva.df$Time.point=="T2",]$Spectral.Slope.Ratio
DOC_T3$sr=suva.df[suva.df$Time.point=="T3",]$Spectral.Slope.Ratio
DOC_T4$sr=suva.df[suva.df$Time.point=="T4",]$Spectral.Slope.Ratio
#freshness (BIX)
DOC_T0$fr=index.df[index.df$Time.point=="T0",]$Freshness.Index
DOC_T1$fr=index.df[index.df$Time.point=="T1",]$Freshness.Index
DOC_T2$fr=index.df[index.df$Time.point=="T2",]$Freshness.Index
DOC_T3$fr=index.df[index.df$Time.point=="T3",]$Freshness.Index
DOC_T4$fr=index.df[index.df$Time.point=="T4",]$Freshness.Index
#combine
DOC.df<-rbind(DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4)
#Extract file
DOC.df$File <- sapply(strsplit(DOC.df$File, "/"), `[`, 3) # extract sample names
# alternative way to code the above
#give the 10th-24th character of the file name, removing the rest
#DOC.df$File<-substr(DOC.df$File, 13, 27)
#alternatively
# remove the 9 letters ('^.) at start
# remove the 4 letters (.$') at end
#DOC.df$File<-gsub('^.........|....$', '', DOC.df$File)
# if equals DOC_T0 then, T0, if not then T1..etc
DOC.df$Time.point<- as.factor(ifelse(DOC.df$File=="DOC_T0.csv", "T0",
ifelse (DOC.df$File=="DOC_T1.csv", "T1",
ifelse (DOC.df$File=="DOC_T2.csv", "T2",
ifelse(DOC.df$File=="DOC_T3.csv", "T3", "T4")))))
#rearrange
library(dplyr)
DOC.df<- DOC.df %>%
dplyr::select(File, Time.point, Treatment, Tank, plant.mass..g, DOC..mg.L, TN..mg.L, SUVA, sr, hix, fl, fr)
DOC.df=DOC.df%>%
mutate_if(is.character,as.factor)SUVA analysis
We can model how SUVA changes across the biomass gradient and over time.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for SUVA.
######## T0 model
m1.SUVA.T0.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.SUVA.T0.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.SUVA.T0.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.suva.AIC<-AIC(m1.SUVA.T0.gam, m2.SUVA.T0.gam, m3.SUVA.T0.gam)
T0.suva.AIC
summary(m1.SUVA.T0.gam)
anova.gam(m1.SUVA.T0.gam)
gam.check(m1.SUVA.T0.gam, rep=1000)
draw(m1.SUVA.T0.gam)
concrvity(m1.SUVA.T0.gam)
par(mfrow = c(2, 2))
plot(m1.SUVA.T0.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T0<-plot_difference(
m1.SUVA.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
suva.T0.mod.plot<-
plot_smooths(
model = m3.SUVA.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
coord_cartesian(ylim = c(0,12))+
xlab("plant material (g)") +
Fig.formatting
suva.T0.mod.plot
######## T1 model
m1.suva.T1.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.suva.T1.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.suva.T1.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.suva.AIC<-AIC(m1.suva.T1.gam, m2.suva.T1.gam, m3.suva.T1.gam)
T1.suva.AIC
summary(m1.suva.T1.gam)
T1.suva.anova=anova(m1.suva.T1.gam)
gam.check(m1.suva.T1.gam, rep=1000)
draw(m1.suva.T1.gam)
concrvity(m1.suva.T1.gam)
par(mfrow = c(2, 2))
plot(m1.suva.T1.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T1<-plot_difference(
m1.suva.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T1.mod.plot<-
plot_smooths(
model = m1.suva.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab(expression('SUVA'["254"]~'(mgC'~ L^-1~m^-1*')'))+
coord_cartesian(ylim = c(0,12))+
xlab("") +
Fig.formatting
########## T2
m1.suva.T2.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.suva.T2.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.suva.T2.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.suva.AIC<-AIC(m1.suva.T2.gam, m2.suva.T2.gam, m3.suva.T2.gam)
T2.suva.AIC
summary(m1.suva.T2.gam)
T2.suva.anova=anova(m1.suva.T2.gam)
gam.check(m1.suva.T2.gam, rep=1000)
draw(m1.suva.T2.gam)
concrvity(m1.suva.T2.gam)
par(mfrow = c(2, 2))
plot(m1.suva.T2.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T2<-plot_difference(
m1.suva.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T2.mod.plot<-
plot_smooths(
model = m1.suva.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0, 12) +
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
suva.T2.mod.plot
########## T3
m1.suva.T3.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.suva.T3.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.suva.T3.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.suva.AIC<-AIC(m1.suva.T3.gam, m2.suva.T3.gam, m3.suva.T3.gam)
T3.suva.AIC
summary(m1.suva.T3.gam)
T3.suva.anova=anova(m1.suva.T3.gam)
gam.check(m1.suva.T3.gam, rep=1000)
draw(m1.suva.T3.gam)
concrvity(m1.suva.T3.gam)
par(mfrow = c(2, 2))
plot(m1.suva.T3.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T3<-plot_difference(
m1.suva.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T3.mod.plot<-
plot_smooths(
model = m1.suva.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,12))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.suva.T4.gam=gam(SUVA ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.suva.T4.gam=gam(SUVA ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.suva.T4.gam=gam(SUVA ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.suva.AIC<-AIC(m1.suva.T4.gam, m2.suva.T4.gam, m3.suva.T4.gam)
T4.suva.AIC
# best is global
summary(m1.suva.T4.gam)
T4.suva.anova=anova(m1.suva.T4.gam)
gam.check(m1.suva.T4.gam, rep=1000)
draw(m1.suva.T4.gam)
concrvity(m1.suva.T4.gam)
par(mfrow = c(2, 2))
plot(m1.suva.T4.gam, all.terms = TRUE, page=1)
# model predictions
suva.diff.T4<-plot_difference(
m1.suva.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-5,5))+
Fig.formatting
###########
#plot for the model output on rawdata
suva.T4.mod.plot<-
plot_smooths(
model = m1.suva.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=SUVA, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,12))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
suva.T4.mod.plot
mod.suva.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.suva.column<- data.frame(mod.suva.rep)
AIC.suva<-bind_rows(T1.suva.AIC, T2.suva.AIC, T3.suva.AIC, T4.suva.AIC)
AIC.suva.mod<-as.data.frame(cbind(mod.suva.column, AIC.suva))
names(AIC.suva.mod)[1]="Model"
AIC.suva.mod <- AIC.suva.mod %>%
# Creating an empty column:
add_column(Variable= c("SUVA"), .before="Model")
write.csv(AIC.suva.mod, "output/AIC models/AIC.suva.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/suva.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.suva.mod, rows= NULL)
dev.off()SUVA output
This generates a table of GAM summaries.
####Parametric output
suva.para.sums=as.data.frame(rbind(T1.suva.anova$pTerms.table, T2.suva.anova$pTerms.table, T3.suva.anova$pTerms.table, T4.suva.anova$pTerms.table))
names(suva.para.sums)[1]="df/edf"
suva.para.sums
#create empty column for ref.df
suva.para.sums <- suva.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
suva.smooth.sums=as.data.frame(rbind(T1.suva.anova$s.table, T2.suva.anova$s.table, T3.suva.anova$s.table, T4.suva.anova$s.table))
names(suva.smooth.sums)[1]="df/edf"
#create empty column for ref.df
suva.smooth.sums <- suva.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
suva.mod.allsums=rbind(suva.para.sums,suva.smooth.sums)
rownames(suva.mod.allsums)<-c(1:12)
#create new ID column
suva.mod.allsums <- suva.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("SUVA"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/suva.mod.allsums.pdf", width = 8.5, height = 3.8) # Export PDF
grid.table(suva.mod.allsums, rows= NULL)
dev.off()SUVA model plot
extract.legend <- get_legend(
# create some space to the left of the legend
DOC.T1.mod.plot + theme(legend.box.margin = margin(0, 0, 0, 10)))
suva.alltimes<-plot_grid(
suva.T1.mod.plot+ theme(legend.position = "none"),
suva.T2.mod.plot+ theme(legend.position = "none"),
suva.T3.mod.plot+ theme(legend.position = "none"),
suva.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
suva.alltimes
Figure. Changes in SUVA with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/SUVA.alltimes.mod.pdf", height=4, width=12)SUVA diff plot
suva.mod.diffs<-plot_grid(
suva.diff.T1+ theme(legend.position = "none")+ ggtitle("SUVA - Day-10"),
suva.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
suva.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
suva.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
suva.mod.diffs
Figure. Differences between fitted smooth functions for the fire treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/suva.mod.diffs.pdf", height=4, width=12)Sr analysis
We can model how Sr changes across the biomass gradient and over time.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for Sr.
######## T0 model
m1.sr.T0.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.sr.T0.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.sr.T0.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.sr.AIC<-AIC(m1.sr.T0.gam, m2.sr.T0.gam, m3.sr.T0.gam)
T0.sr.AIC
summary(m1.sr.T0.gam)
anova.gam(m1.sr.T0.gam)
gam.check(m1.sr.T0.gam, rep=1000)
draw(m1.sr.T0.gam)
concrvity(m1.sr.T0.gam)
par(mfrow = c(2, 2))
plot(m1.sr.T0.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T0<-plot_difference(
m1.sr.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
sr.T0.mod.plot<-
plot_smooths(
model = m3.sr.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylab(expression('S'[italic(R)]))+
xlab("plant material (g)") +
Fig.formatting
sr.T0.mod.plot
######## T1 model
m1.sr.T1.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.sr.T1.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.sr.T1.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.sr.AIC<-AIC(m1.sr.T1.gam, m2.sr.T1.gam, m3.sr.T1.gam)
T1.sr.AIC
summary(m1.sr.T1.gam)
T1.sr.anova=anova(m1.sr.T1.gam)
gam.check(m1.sr.T1.gam, rep=1000)
draw(m1.sr.T1.gam)
concrvity(m1.sr.T1.gam)
par(mfrow = c(2, 2))
plot(m1.sr.T1.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T1<-plot_difference(
m1.sr.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T1.mod.plot<-
plot_smooths(
model = m1.sr.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab(expression('S'[italic(R)]))+
coord_cartesian(ylim = c(0.25,2.5))+
xlab("") +
Fig.formatting
sr.T1.mod.plot
########## T2
m1.sr.T2.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.sr.T2.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.sr.T2.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.sr.AIC<-AIC(m1.sr.T2.gam, m2.sr.T2.gam, m3.sr.T2.gam)
T2.sr.AIC
summary(m3.sr.T2.gam)
T2.sr.anova=anova(m1.sr.T2.gam)
gam.check(m3.sr.T2.gam, rep=1000)
draw(m3.sr.T2.gam)
concrvity(m3.sr.T2.gam)
par(mfrow = c(2, 2))
plot(m3.sr.T2.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T2<-plot_difference(
m1.sr.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T2.mod.plot<-
plot_smooths(
model = m3.sr.T2.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.25,2.5))+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
sr.T2.mod.plot
########## T3
m1.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.sr.T3.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m4.sr.T3.gam=gam(sr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "ML")
T3.sr.AIC<-AIC(m1.sr.T3.gam, m2.sr.T3.gam, m3.sr.T3.gam)
T3.sr.AIC
summary(m1.sr.T3.gam)
T3.sr.anova=anova(m1.sr.T3.gam)
gam.check(m1.sr.T3.gam, rep=1000)
draw(m1.sr.T3.gam)
concrvity(m1.sr.T3.gam)
par(mfrow = c(2, 2))
plot(m1.sr.T3.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T3<-plot_difference(
m1.sr.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T3.mod.plot<-
plot_smooths(
model = m2.sr.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen"))+
coord_cartesian(ylim = c(0.25,2.5))+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
sr.T3.mod.plot
########## T4
m1.sr.T4.gam=gam(sr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.sr.T4.gam=gam(sr ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.sr.T4.gam=gam(sr ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.sr.AIC<-AIC(m1.sr.T4.gam, m2.sr.T4.gam, m3.sr.T4.gam)
T4.sr.AIC
summary(m1.sr.T4.gam)
T4.sr.anova=anova(m1.sr.T4.gam)
gam.check(m1.sr.T4.gam, rep=1000)
draw(m1.sr.T4.gam)
concrvity(m1.sr.T4.gam)
par(mfrow = c(2, 2))
plot(m1.sr.T4.gam, all.terms = TRUE, page=1)
# model predictions
sr.diff.T4<-plot_difference(
m1.sr.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-0.75,0.75))+
Fig.formatting
###########
#plot for the model output on rawdata
sr.T4.mod.plot<-
plot_smooths(
model = m1.sr.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=sr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.25,2.5))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
sr.T4.mod.plot
mod.sr.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.sr.column<- data.frame(mod.sr.rep)
AIC.sr<-bind_rows(T1.sr.AIC, T2.sr.AIC, T3.sr.AIC, T4.sr.AIC)
AIC.sr
AIC.sr.mod<-as.data.frame(cbind(mod.sr.column, AIC.sr))
names(AIC.sr.mod)[1]="Model"
AIC.sr.mod <- AIC.sr.mod %>%
# Creating an empty column:
add_column(Variable= c("Slope ratio"), .before="Model")
write.csv(AIC.sr.mod, "output/AIC models/AIC.sr.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/sr.mod.AIC.pdf", width = 6.3, height = 3.8) # Export PDF
grid.table(AIC.sr.mod, rows=NULL)
dev.off()Sr output
This generates a table GAM summaries.
####Parametric output
sr.para.sums=as.data.frame(rbind(T1.sr.anova$pTerms.table, T2.sr.anova$pTerms.table, T3.sr.anova$pTerms.table, T4.sr.anova$pTerms.table))
names(sr.para.sums)[1]="df/edf"
#create empty column for ref.df
sr.para.sums <- sr.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
sr.smooth.sums=as.data.frame(rbind(T1.sr.anova$s.table, T2.sr.anova$s.table, T3.sr.anova$s.table, T4.sr.anova$s.table))
names(sr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
sr.smooth.sums <- sr.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
sr.smooth.sums
#merge the dataframes
sr.mod.allsums=rbind(sr.para.sums,sr.smooth.sums)
rownames(sr.mod.allsums)<-c(1:12)
#create new ID column
sr.mod.allsums <- sr.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("Slope ratio"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/sr.mod.allsums.pdf", width = 8.5, height = 3.8) # Export PDF
grid.table(sr.mod.allsums, rows= NULL)
dev.off()Sr model plot
sr.alltimes<-plot_grid(
sr.T1.mod.plot+ theme(legend.position = "none"),
sr.T2.mod.plot+ theme(legend.position = "none"),
sr.T3.mod.plot+ theme(legend.position = "none"),
sr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
sr.alltimes
Figure. Changes in Sr with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/sr.alltimes.mod.pdf", height=4, width=12)Sr diff plot
sr.mod.diffs<-plot_grid(
sr.diff.T1+ theme(legend.position = "none")+ ggtitle("Sr - Day-10"),
sr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
sr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
sr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
sr.mod.diffs
Figure. Differences between fitted smooth functions for the fire treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/sr.mod.diffs.pdf", height=4, width=12)HIX analysis
We can model how HIX changes across the biomass gradient and over time.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for HIX.
######## T0 model
m1.hix.T0.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.hix.T0.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.hix.T0.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.hix.AIC<-AIC(m1.hix.T0.gam, m2.hix.T0.gam, m3.hix.T0.gam)
T0.hix.AIC
summary(m1.hix.T0.gam)
anova.gam(m1.hix.T0.gam)
gam.check(m1.hix.T0.gam, rep=1000)
draw(m1.hix.T0.gam)
concrvity(m1.hix.T0.gam)
par(mfrow = c(2, 2))
plot(m1.hix.T0.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T0<-plot_difference(
m1.hix.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
hix.T0.mod.plot<-
plot_smooths(
model = m3.hix.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
coord_cartesian(ylim = c(0,6))+
ylab("HIX")+
xlab("plant material (g)") +
Fig.formatting
######## T1 model
m1.hix.T1.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.hix.T1.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.hix.T1.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.hix.AIC<-AIC(m1.hix.T1.gam, m2.hix.T1.gam, m3.hix.T1.gam)
T1.hix.AIC
summary(m3.hix.T1.gam)
T1.hix.anova=anova(m1.hix.T1.gam)
gam.check(m3.sr.T1.gam, rep=1000)
draw(m3.hix.T1.gam)
concrvity(m3.hix.T1.gam)
par(mfrow = c(2, 2))
plot(m3.hix.T1.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T1<-plot_difference(
m1.hix.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T1.mod.plot<-
plot_smooths(
model = m3.hix.T1.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-10") +
ylab("HIX")+
xlab("") +
Fig.formatting
hix.T1.mod.plot
########## T2
m1.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.hix.T2.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m4.hix.T2.gam=gam(hix ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.hix.AIC<-AIC(m1.hix.T2.gam, m2.hix.T2.gam, m3.hix.T2.gam)
T2.hix.AIC
summary(m2.hix.T2.gam)
T2.hix.anova=anova(m1.hix.T2.gam)
gam.check(m2.hix.T2.gam, rep=1000)
draw(m2.hix.T2.gam)
concrvity(m2.hix.T2.gam)
par(mfrow = c(2, 2))
plot(m2.hix.T2.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T2<-plot_difference(
m1.hix.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T2.mod.plot<-
plot_smooths(
model = m2.hix.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
hix.T2.mod.plot
########## T3
m1.hix.T3.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.hix.T3.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.hix.T3.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.hix.AIC<-AIC(m1.hix.T3.gam, m2.hix.T3.gam, m3.hix.T3.gam)
T3.hix.AIC
summary(m1.hix.T3.gam)
T3.hix.anova=anova(m1.hix.T3.gam)
gam.check(m1.hix.T3.gam, rep=1000)
draw(m1.hix.T3.gam)
concrvity(m1.hix.T3.gam)
par(mfrow = c(2, 2))
plot(m1.hix.T3.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T3<-plot_difference(
m1.hix.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T3.mod.plot<-
plot_smooths(
model = m1.hix.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.hix.T4.gam=gam(hix ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.hix.T4.gam=gam(hix ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.hix.T4.gam=gam(hix ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.hix.AIC<-AIC(m1.hix.T4.gam, m2.hix.T4.gam, m3.hix.T4.gam)
T4.hix.AIC
summary(m1.hix.T4.gam)
T4.hix.anova=anova(m1.hix.T4.gam)
gam.check(m1.hix.T4.gam, rep=1000)
draw(m1.hix.T4.gam)
concrvity(m1.hix.T4.gam)
par(mfrow = c(2, 2))
plot(m1.hix.T4.gam, all.terms = TRUE, page=1)
# model predictions
hix.diff.T4<-plot_difference(
m1.hix.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-2.5,1.5))+
Fig.formatting
###########
#plot for the model output on rawdata
hix.T4.mod.plot<-
plot_smooths(
model = m1.hix.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=hix, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0,6))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
mod.hix.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.hix.column<- data.frame(mod.hix.rep)
AIC.hix<-bind_rows(T1.hix.AIC, T2.hix.AIC, T3.hix.AIC, T4.hix.AIC)
AIC.hix.mod<-as.data.frame(cbind(mod.hix.column, AIC.hix))
names(AIC.hix.mod)[1]="Model"
AIC.hix.mod <- AIC.hix.mod %>%
# Creating an empty column:
add_column(Variable= c("HIX"), .before="Model")
write.csv(AIC.hix.mod, "output/AIC models/AIC.hix.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/hix.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.hix.mod, rows= NULL)
dev.off()HIX output
This generates a table with GAM summaries.
####Parametric output
hix.para.sums=as.data.frame(rbind(T1.hix.anova$pTerms.table, T2.hix.anova$pTerms.table, T3.hix.anova$pTerms.table, T4.hix.anova$pTerms.table))
names(hix.para.sums)[1]="df/edf"
#create empty column for ref.df
hix.para.sums <- hix.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
hix.smooth.sums=as.data.frame(rbind(T1.hix.anova$s.table, T2.hix.anova$s.table, T3.hix.anova$s.table, T4.hix.anova$s.table))
names(hix.smooth.sums)[1]="df/edf"
#create empty column for ref.df
hix.smooth.sums <- hix.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
hix.mod.allsums=rbind(hix.para.sums,hix.smooth.sums)
#create new ID column
hix.mod.allsums <- hix.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("HIX"), .before="Model")
#convert dataframe into a table for export
pdf("output/Codys.tables/hix.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(hix.mod.allsums, rows= NULL)
dev.off()HIX model plot
hix.alltimes<-plot_grid(
hix.T1.mod.plot+ theme(legend.position = "none"),
hix.T2.mod.plot+ theme(legend.position = "none"),
hix.T3.mod.plot+ theme(legend.position = "none"),
hix.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
hix.alltimes
Figure. Changes in HIX with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/hix.alltimes.mod.pdf", height=4, width=12)HIX diff plot
hix.mod.diffs<-plot_grid(
hix.diff.T1+ theme(legend.position = "none")+ ggtitle("HIX - Day-10"),
hix.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
hix.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
hix.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
hix.mod.diffs
Figure. Differences between fitted smooth functions for the fire treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/hix.mod.diffs.pdf", height=4, width=12)BIX analysis
We can model how BIX changes across the biomass gradient and over time.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for BIX.
######## T0 model
m1.fr.T0.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.fr.T0.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.fr.T0.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.fr.AIC<-AIC(m1.fr.T0.gam, m2.fr.T0.gam, m3.fr.T0.gam)
T0.fr.AIC
summary(m1.fr.T0.gam)
anova.gam(m1.fr.T0.gam)
gam.check(m1.fr.T0.gam, rep=1000)
draw(m1.fr.T0.gam)
concrvity(m1.fr.T0.gam)
par(mfrow = c(2, 2))
plot(m1.fr.T0.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T0<-plot_difference(
m1.fr.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)
###########
fr.T0.mod.plot<-
plot_smooths(
model = m3.fr.T0.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
ylim(0.25,1.75)+
xlab("plant material (g)") +
Fig.formatting
fr.T0.mod.plot
######## T1 model
m1.fr.T1.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.fr.T1.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.fr.T1.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.fr.AIC<-AIC(m1.fr.T1.gam, m2.fr.T1.gam, m3.fr.T1.gam)
T1.fr.AIC
summary(m1.fr.T1.gam)
T1.fr.anova=anova(m1.fr.T1.gam)
gam.check(m1.fr.T1.gam, rep=1000)
draw(m1.fr.T1.gam)
concrvity(m1.fr.T1.gam)
par(mfrow = c(2, 2))
plot(m1.fr.T1.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T1<-plot_difference(
m1.fr.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T1.mod.plot<-
plot_smooths(
model = m1.fr.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylim(0.25,1.75)+
ylab(expression("Freshness"~"("~italic(alpha)~":"~italic(beta)~")"))+
xlab("") +
Fig.formatting
fr.T1.mod.plot
########## T2
m1.fr.T2.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.fr.T2.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.fr.T2.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.fr.AIC<-AIC(m1.fr.T2.gam, m2.fr.T2.gam, m3.fr.T2.gam)
T2.fr.AIC
summary(m1.fr.T2.gam)
T2.fr.anova=anova(m1.fr.T2.gam)
gam.check(m1.fr.T2.gam, rep=1000)
draw(m1.fr.T2.gam)
concrvity(m1.fr.T2.gam)
par(mfrow = c(2, 2))
plot(m1.fr.T2.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T2<-plot_difference(
m1.fr.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T2.mod.plot<-
plot_smooths(
model = m1.fr.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
Fig.formatting
fr.T2.mod.plot
########## T3
m1.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.fr.T3.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m4.fr.T3.gam=gam(fr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.fr.AIC<-AIC(m1.fr.T3.gam, m2.fr.T3.gam, m3.fr.T3.gam)
T3.fr.AIC
summary(m2.fr.T3.gam)
T3.fr.anova=anova(m1.fr.T3.gam)
gam.check(m2.fr.T3.gam, rep=1000)
draw(m2.fr.T3.gam)
concrvity(m2.fr.T3.gam)
par(mfrow = c(2, 2))
plot(m2.fr.T3.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T3<-plot_difference(
m1.fr.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T3.mod.plot<-
plot_smooths(
model = m2.fr.T3.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
geom_line(aes(linetype=Treatment))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.fr.T4.gam=gam(fr ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m4.fr.T4.gam=gam(fr ~ Treatment + s(plant.mass..g, m=c(2,0)), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.fr.AIC<-AIC(m1.fr.T4.gam, m2.fr.T4.gam, m3.fr.T4.gam)
T4.fr.AIC
summary(m2.fr.T4.gam)
T4.fr.anova=anova(m1.fr.T4.gam)
gam.check(m2.fr.T4.gam, rep=1000)
draw(m2.fr.T4.gam)
concrvity(m2.fr.T4.gam)
par(mfrow = c(2, 2))
plot(m2.fr.T4.gam, all.terms = TRUE, page=1)
# model predictions
fr.diff.T4<-plot_difference(
m1.fr.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1,0.2))+
Fig.formatting
###########
#plot for the model output on rawdata
fr.T4.mod.plot<-
plot_smooths(
model = m2.fr.T4.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=fr, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ylim(0.25,1.75)+
ggtitle("Day-89") +
geom_line(aes(linetype=Treatment))+
ylab("") +
xlab("") +
Fig.formattingBIX output
This generates a table with GAM summaries.
####Parametric output
fr.para.sums=as.data.frame(rbind(T1.fr.anova$pTerms.table, T2.fr.anova$pTerms.table, T3.fr.anova$pTerms.table, T4.fr.anova$pTerms.table))
names(fr.para.sums)[1]="df/edf"
#create empty column for ref.df
fr.para.sums <- fr.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
fr.smooth.sums=as.data.frame(rbind(T1.fr.anova$s.table, T2.fr.anova$s.table, T3.fr.anova$s.table, T4.fr.anova$s.table))
names(fr.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fr.smooth.sums <- fr.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
fr.mod.allsums=rbind(fr.para.sums,fr.smooth.sums)
#create new ID column
fr.mod.allsums <- fr.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("FR"), .before="Model")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(fr.mod.allsums, rows= NULL)
dev.off()
mod.fr.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.fr.column<- data.frame(mod.fr.rep)
AIC.fr<-bind_rows(T1.fr.AIC, T2.fr.AIC, T3.fr.AIC, T4.fr.AIC)
AIC.fr.mod<-as.data.frame(cbind(mod.fr.column, AIC.fr))
names(AIC.fr.mod)[1]="Model"
AIC.fr.mod <- AIC.fr.mod %>%
# Creating an empty column:
add_column(Variable= c("FR"), .before="Model")
write.csv(AIC.fr.mod, "output/AIC models/AIC.fr.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fr.mod.AIC.pdf", width = 6.3, height = 3.8) # Export PDF
grid.table(AIC.fr.mod, rows= NULL)
dev.off()BIX model plot
fr.alltimes<-plot_grid(
fr.T1.mod.plot+ theme(legend.position = "none"),
fr.T2.mod.plot+ theme(legend.position = "none"),
fr.T3.mod.plot+ theme(legend.position = "none"),
fr.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
fr.alltimes
Figure. Changes in BIX with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)BIX diff plot
fr.mod.diffs<-plot_grid(
fr.diff.T1+ theme(legend.position = "none")+ ggtitle("FR - Day-10"),
fr.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
fr.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
fr.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
fr.mod.diffs
Figure. Differences between fitted smooth functions for the fire treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/fr.mod.diffs.pdf", height=4, width=12)Fl analysis
We can model how Fl changes across the biomass gradient and over time.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for Fl.
library(mgcv)
library(tidymv)
######## T0 model
m1.fl.T0.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T0", data = DOC.df, method = "REML")
m2.fl.T0.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
m3.fl.T0.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T0", data = DOC.df, method = "REML")
T0.fl.AIC<-AIC(m1.fl.T0.gam, m2.fl.T0.gam, m3.fl.T0.gam)
T0.fl.AIC
summary(m1.fl.T0.gam)
anova.gam(m1.fl.T0.gam)
gam.check(m1.fl.T0.gam, rep=1000)
draw(m1.fl.T0.gam)
concrvity(m1.fl.T0.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T0.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T0<-plot_difference(
m1.fl.T0.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,0.5))+
Fig.formatting
###########
fl.T0.mod.plot<-
plot_smooths(
model = m1.fl.T0.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T0"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-0") +
coord_cartesian(ylim = c(0.5,2.5))+
ylab("FI")+
xlab("plant material (g)") +
Fig.formatting
fl.T0.mod.plot
######## T1 model
m1.fl.T1.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T1", data = DOC.df, method = "REML")
m2.fl.T1.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
m3.fl.T1.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T1", data = DOC.df, method = "REML")
T1.fl.AIC<-AIC(m1.fl.T1.gam, m2.fl.T1.gam, m3.fl.T1.gam)
T1.fl.AIC
summary(m1.fl.T1.gam)
T1.fl.anova=anova(m1.fl.T1.gam)
T1.fl.anova
gam.check(m1.fl.T1.gam, rep=1000)
draw(m1.fl.T1.gam)
concrvity(m1.fl.T1.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T1.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T1<-plot_difference(
m1.fl.T1.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T1.mod.plot<-
plot_smooths(
model = m1.fl.T1.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-10") +
ylab("FI")+
xlab("") +
ylim(0.5,2.5)+
Fig.formatting
fl.T1.mod.plot
########## T2
m1.fl.T2.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T2", data = DOC.df, method = "REML")
m2.fl.T2.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
m3.fl.T2.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T2", data = DOC.df, method = "REML")
T2.fl.AIC<-AIC(m1.fl.T2.gam, m2.fl.T2.gam, m3.fl.T2.gam)
T2.fl.AIC
summary(m1.fl.T2.gam)
T2.fl.anova=anova(m1.fl.T2.gam)
gam.check(m1.fl.T2.gam, rep=1000)
draw(m1.fl.T2.gam)
concrvity(m1.fl.T2.gam)
par(mfrow = c(2, 2))
plot(m1.fl.T2.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T2<-plot_difference(
m1.fl.T2.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T2.mod.plot<-
plot_smooths(
model = m1.fl.T2.gam,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T2"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Day-31") +
ylab("") +
xlab("plant material (g)") +
coord_cartesian(ylim = c(0.5,2.5))+
Fig.formatting
fl.T2.mod.plot
########## T3
m1.fl.T3.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T3", data = DOC.df, method = "REML")
m2.fl.T3.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
m3.fl.T3.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T3", data = DOC.df, method = "REML")
T3.fl.AIC<-AIC(m1.fl.T3.gam, m2.fl.T3.gam, m3.fl.T3.gam)
T3.fl.AIC
summary(m3.fl.T3.gam)
T3.fl.anova=anova(m1.fl.T3.gam)
gam.check(m3.fl.T3.gam, rep=1000)
draw(m3.fl.T3.gam)
concrvity(m3.fl.T3.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T3.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T3<-plot_difference(
m1.fl.T3.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T3.mod.plot<-
plot_smooths(
model = m3.fl.T3.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T3"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
ggtitle("Day-59") +
ylab("") +
xlab("") +
Fig.formatting
########## T4
m1.fl.T4.gam=gam(fl ~ Treatment + s(plant.mass..g, by= Treatment), subset = Time.point=="T4", data = DOC.df, method = "REML")
m2.fl.T4.gam=gam(fl ~ Treatment + s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
m3.fl.T4.gam=gam(fl ~ s(plant.mass..g), subset = Time.point=="T4", data = DOC.df, method = "REML")
T4.fl.AIC<-AIC(m1.fl.T4.gam, m2.fl.T4.gam, m3.fl.T4.gam)
T4.fl.AIC
summary(m3.fl.T4.gam)
T4.fl.anova=anova(m1.fl.T4.gam)
gam.check(m3.fl.T4.gam, rep=1000)
draw(m3.fl.T4.gam)
concrvity(m3.fl.T4.gam)
par(mfrow = c(2, 2))
plot(m3.fl.T4.gam, all.terms = TRUE, page=1)
# model predictions
fl.diff.T4<-plot_difference(
m1.fl.T4.gam,
series = plant.mass..g,
difference = list(Treatment = c("burned", "unburned"))
)+
coord_cartesian(ylim = c(-1.5,1))+
Fig.formatting
###########
#plot for the model output on rawdata
fl.T4.mod.plot<-
plot_smooths(
model = m3.fl.T4.gam,
series = plant.mass..g,
) +
geom_point(data=DOC.df[(DOC.df$Time.point=="T4"),],
aes(x=plant.mass..g, y=fl, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
coord_cartesian(ylim = c(0.5,2.5))+
ggtitle("Day-89") +
ylab("") +
xlab("") +
Fig.formatting
mod.fl.rep<-rep(c(
"Day 10: Treatment with by-factor smooth",
"Day 10: Treatment with global smooth",
"Day 10: global smooth",
"Day 31: Treatment with by-factor smooth",
"Day 31: Treatment with global smooth",
"Day 31: global smooth",
"Day 59: Treatment with by-factor smooth",
"Day 59: Treatment with global smooth",
"Day 59: global smooth",
"Day 89: Treatment with by-factor smooth",
"Day 89: Treatment with global smooth",
"Day 89: global smooth"))
mod.fl.column<- data.frame(mod.fl.rep)
AIC.fl<-bind_rows(T1.fl.AIC, T2.fl.AIC, T3.fl.AIC, T4.fl.AIC)
AIC.fl.mod<-as.data.frame(cbind(mod.fl.column, AIC.fl))
AIC.fl.mod
names(AIC.fl.mod)[1]="Model"
AIC.fl.mod <- AIC.fl.mod %>%
# Creating an empty column:
add_column(Variable= c("FI"), .before="Model")
write.csv(AIC.fl.mod, "output/AIC models/AIC.fl.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.AIC.pdf", width = 6, height = 3.8) # Export PDF
grid.table(AIC.fl.mod, rows=NULL)
dev.off()Fl output
This generates a table with GAM summaries.
####Parametric output
fl.para.sums=as.data.frame(rbind(T1.fl.anova$pTerms.table, T2.fl.anova$pTerms.table, T3.fl.anova$pTerms.table, T4.fl.anova$pTerms.table))
names(fl.para.sums)[1]="df/edf"
#create empty column for ref.df
fl.para.sums <- fl.para.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: Treatment","Day 31 Model: Treatment","Day 59 Model: Treatment","Day 89 Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
fl.smooth.sums=as.data.frame(rbind(T1.fl.anova$s.table, T2.fl.anova$s.table, T3.fl.anova$s.table, T4.fl.anova$s.table))
names(fl.smooth.sums)[1]="df/edf"
#create empty column for ref.df
fl.smooth.sums <- fl.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Day 10 Model: burned smooth","Day 10 Model: unburned smooth","Day 31 Model: burned smooth","Day 31 Model: unburned smooth","Day 59 Model: burned smooth","Day 59 Model: unburned smooth","Day 89 Model: burned smooth","Day 89 Model: unburned smooth"), .before="df/edf")
#merge the dataframes
fl.mod.allsums=rbind(fl.para.sums,fl.smooth.sums)
rownames(fl.mod.allsums)<-c(1:12)
#create new ID column
fl.mod.allsums <- fl.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("FI"), .before="Model")
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/fl.mod.allsums.pdf", width = 8.3, height = 3.8) # Export PDF
grid.table(fl.mod.allsums, rows=NULL)
dev.off()Fl model plot
fl.alltimes<-plot_grid(
fl.T1.mod.plot+ theme(legend.position = "none"),
fl.T2.mod.plot+ theme(legend.position = "none"),
fl.T3.mod.plot+ theme(legend.position = "none"),
fl.T4.mod.plot+ theme(legend.position = "none"), extract.legend,
rel_widths = c(8,8,8,8,3), ncol=5)
fl.alltimes
Figure. Changes in Fl with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/fr.alltimes.mod.pdf", height=4, width=12)Fl diff plot
fl.mod.diffs<-plot_grid(
fl.diff.T1+ theme(legend.position = "none")+ ggtitle("FI - Day-10"),
fl.diff.T2+ theme(legend.position = "none")+ ggtitle("Day-31"),
fl.diff.T3+ theme(legend.position = "none")+ ggtitle("Day-59"),
fl.diff.T4+ theme(legend.position = "none")+ ggtitle("Day-89"),
rel_widths = c(8,8,8,8), ncol=4)
fl.mod.diffs
Figure. Differences between fitted smooth functions for the fire treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/fl.mod.diffs.pdf", height=4, width=12)All indices model plot
This is a large 5 x 4 figure with all fluorescence and UV-vis absorbance indices across sampling time points.
indices.alltimes= plot_grid(suva.alltimes, sr.alltimes, hix.alltimes, fr.alltimes, fl.alltimes, nrow=5, rel_widths = c(8,8,8,8,8), labels = "AUTO")
indices.alltimes
Figure. Changes in a: SUVA, b: spectral slope ratio, c: humification index, d: freshness index, and e: fluorescence index with increasing plant biomass for the burned (red) and unburned (green) treatments. Each dot represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/indices.alltimes.pdf", height=18, width=14)
#### Day 0 baseline indicies
day0_index.baseline.plot = plot_grid(
suva.T0.mod.plot+ theme(legend.position = "none"),
sr.T0.mod.plot+ theme(legend.position = "none"),
hix.T0.mod.plot+ theme(legend.position = "none"),
fr.T0.mod.plot+ theme(legend.position = "none"),
fl.T0.mod.plot+ theme(legend.position = "none"),
extract.legend,
rel_widths = c(8,8,8,8,8,3), nrow=6)
ggsave("figures/day0_index.baseline.plot.pdf", height=18, width=14)All indices diff plot
indices.diff.plot= plot_grid(suva.mod.diffs, sr.mod.diffs, hix.mod.diffs, fr.mod.diffs, fl.mod.diffs, nrow=5, rel_widths = c(8,8,8,8,8), labels = "auto")
indices.diff.plot
Figure.Differences between fitted smooth functions for the treatments (difference in trends; solid lines) for fluorescence and UV-vis absorbance indices shown in red with approximate 95% confidence intervals.
ggsave("figures/indices.diff.plot.pdf", height=18, width=14)Import UVBIO data
Import the data for UVBIO (data from the in situ incubations) and run a loop to clean up all files. This creates a single df by taking the raw data files, aligning metadata, and filtering to make a clean df for GAMs.
- We will analyze % DOC change from the incubations at Day 10.
- Run model selection and produce plots for each degradation pathway (UV-only, Microbes-only, UV and Microbes).
- Plot the UVBIO data.
##### grab files in a list
UVBio.files <- list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = T)
UVBio.files
##### what are the file names, sans extensions using package 'tools'
file.names<-file_path_sans_ext(list.files(path="/Users/Cody1/Desktop/GitHub/Pyromania.Fiery_death/data/UVBIO", pattern = "csv$", full.names = F))
file.names
############ formatting all data in for loop
for(i in 1:length(UVBio.files))
{
data<-read.csv(UVBio.files[i], sep=",") #skip 13 rows where standards are
data<-data[,c(1,3,4)] # removed columns we don't need
data$File<-UVBio.files[i]
colnames(data)<-c("Tank", "DOC..mg.L", "TN..mg.L", "File")
data$Tank<-as.character(data$Tank)
data$Treatment<-IDs$Treatment # add treatment info
data$plant.mass..g<-IDs$plant.mass..g # add plant biomass info
make.names(assign(paste(file.names[i], sep=""), data)) # make the file name the name of new df for loop df
}
########## this is the end of the loop
ls() #see all dfs you've made, the above will be df matching their file names
#Combine files from loop to single df
UVBIO.df<-rbind(DOC_NoUV_T1, DOC_UV_T1, DOC_NoUV_T2, DOC_UV_T2, DOC_NoUV_T3, DOC_UV_T3)
# make new factors for filtered (F/NF) and UV (+/-)
# ** caution! Make sure this matches sheet layout after the loop!
UVBIO.df$Filter.Trt<-rep(c("Filtered", "Unfiltered"), each=30) # repeat 30 times
UVBIO.df$UV.Trt<-rep(c("NoUV", "UV"), each=60) # repeat 60 times
# check tank names in DF and make sure they are what you think they are (i.e 1-30 but w/ weird names)
Tank.fix<-rep(c(1:30), times=12) # repeat 12 times to fix all the strange tank names
#compare tank IDS
df.test<-cbind(Tank.fix, UVBIO.df$Tank) # yep!
UVBIO.df$Tank<-as.factor(Tank.fix) #replace old tank names
# extract out the file names
#Extract file
UVBIO.df$File=sapply(strsplit(UVBIO.df$File, "/"), `[`, 9)
#take out redundancies
UVBIO.df$DOC..mg.L<-gsub("NPOC:","", UVBIO.df$DOC..mg.L)
UVBIO.df$DOC..mg.L<- gsub("mg/L", "", UVBIO.df$DOC..mg.L)
UVBIO.df$TN..mg.L<-gsub("TN:","", UVBIO.df$TN..mg.L)
UVBIO.df$TN..mg.L<- gsub("mg/L", "", UVBIO.df$TN..mg.L)
#convert DOC and TN to numeric
UVBIO.df$DOC..mg.L<-as.numeric(as.character(UVBIO.df$DOC..mg.L))
UVBIO.df$TN..mg.L<- as.numeric(as.character(UVBIO.df$TN..mg.L))
# add a time point
# if equals DOC_UV_T1_111721 then, T1, if not then "nothing"
UVBIO.df$Time.point<- as.factor(ifelse(UVBIO.df$File=="DOC_UV_T1.csv", "T1",
ifelse(UVBIO.df$File=="DOC_NoUV_T1.csv", "T1",
ifelse(UVBIO.df$File=="DOC_UV_T2.csv", "T2",
ifelse(UVBIO.df$File=="DOC_NoUV_T2.csv", "T2", "T3")))))
#rearrange
UVBIO.df<- UVBIO.df %>%
dplyr::select(File, Time.point, Treatment, Filter.Trt, UV.Trt, Tank, plant.mass..g, DOC..mg.L, TN..mg.L)
# make levels for plotting and stats
# add in starting DOC and TN
# original data is... DOC_T0, DOC_T1, DOC_T2, DOC_T3, DOC_T4 -- but no need for DOC T0 and T4
# repeat 4x to account for the 4 treatments in each tank
DOC_T1.start<-rep(DOC_T1$DOC..mg.L, times=4)
DOC_T2.start<-rep(DOC_T2$DOC..mg.L, times=4)
DOC_T3.start<-rep(DOC_T3$DOC..mg.L, times=4)
TN_T1.start<-rep(DOC_T1$TN..mg.L, times=4)
TN_T2.start<-rep(DOC_T2$TN..mg.L, times=4)
TN_T3.start<-rep(DOC_T3$TN..mg.L, times=4)
# compile and add in as a new column
UVBIO.df$Start.DOC<-c(DOC_T1.start, DOC_T2.start, DOC_T3.start)
UVBIO.df$Start.TN<-c(TN_T1.start, TN_T2.start, TN_T3.start)
UVBIO.df$Burn.Filt.UV<-interaction(UVBIO.df$Treatment, UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)
UVBIO.df$Filt.UV<-interaction(UVBIO.df$Filter.Trt, UVBIO.df$UV.Trt)Create UVBIO df
Since we initially imported the UVBIO files, we can now create a unique df for each degradation pathway and merge them: photo.df = UV-only, microbes.df= microbes-only, total.df= UV and microbes. The degra.df now contains all 3 pathways.
#+UV/+Bio=A
#-UV/+Bio=B
#+UV/-Bio=C
#-UV/-Bio (NoUV,Filtered)=D
#This is calculating % decrease ((control-treatment)/control)*100
#photodegradation
photo.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L),
UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Filtered.UV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
photo.df$effect="UV-only"
names(photo.df)[1]="Time.point"
names(photo.df)[2]="Start.DOC"
names(photo.df)[3]="Start.TN"
names(photo.df)[4]="Treatment"
names(photo.df)[5]="Tank"
names(photo.df)[6]="plant.mass..g"
names(photo.df)[7]="abs.change.DOC"
names(photo.df)[8]="percent.change.DOC"
names(photo.df)[9]="abs.change.TN"
names(photo.df)[10]="percent.change.TN"
names(photo.df)[11]="SUVA"
names(photo.df)[12]="Spectral.slope"
names(photo.df)[13]="HIX"
names(photo.df)[14]="Fr"
names(photo.df)[15]="FI"
photo.df$mass.remaining.DOC=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.DOC)
photo.df$mass.remaining.TN=100-abs(photo.df[photo.df$effect=="UV-only",]$percent.change.TN)
photo_t1.df=subset(photo.df, Time.point %in% c("T1"))
photo.df=photo.df%>%
mutate_if(is.character,as.factor)
photo_t1.df=photo_t1.df%>%
mutate_if(is.character,as.factor)
#Microbe decomp
microbe.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L- UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(( UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)/((UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.NoUV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
microbe.df$effect="microbes-only"
names(microbe.df)[1]="Time.point"
names(microbe.df)[2]="Start.DOC"
names(microbe.df)[3]="Start.TN"
names(microbe.df)[4]="Treatment"
names(microbe.df)[5]="Tank"
names(microbe.df)[6]="plant.mass..g"
names(microbe.df)[7]="abs.change.DOC"
names(microbe.df)[8]="percent.change.DOC"
names(microbe.df)[9]="abs.change.TN"
names(microbe.df)[10]="percent.change.TN"
names(microbe.df)[11]="SUVA"
names(microbe.df)[12]="Spectral.slope"
names(microbe.df)[13]="HIX"
names(microbe.df)[14]="Fr"
names(microbe.df)[15]="FI"
microbe.df$mass.remaining.DOC=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.DOC)
microbe.df$mass.remaining.TN=100-abs(microbe.df[microbe.df$effect=="microbes-only",]$percent.change.TN)
microbe_t1.df=subset(microbe.df, Time.point %in% c("T1"))
microbe.df=microbe.df%>%
mutate_if(is.character,as.factor)
microbe_t1.df=microbe_t1.df%>%
mutate_if(is.character,as.factor)
#Total decomp
total.df=data.frame(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Time.point,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.DOC,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Treatment,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Tank,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$plant.mass..g,
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$DOC..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Filtered.NoUV",]$DOC..mg.L)),
UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN,
100*(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$TN..mg.L-UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)/(abs(UVBIO.df[UVBIO.df$Filt.UV=="Unfiltered.UV",]$Start.TN)),
suva.df[suva.df$Time.point=="T1",]$SUVA,
suva.df[suva.df$Time.point=="T1",]$Spectral.Slope.Ratio,
index.df[index.df$Time.point=="T1",]$Humification.Index,
index.df[index.df$Time.point=="T1",]$Freshness.Index,
index.df[index.df$Time.point=="T1",]$Fluroscence.Index)
total.df$effect="UV and microbes"
names(total.df)[1]="Time.point"
names(total.df)[2]="Start.DOC"
names(total.df)[3]="Start.TN"
names(total.df)[4]="Treatment"
names(total.df)[5]="Tank"
names(total.df)[6]="plant.mass..g"
names(total.df)[7]="abs.change.DOC"
names(total.df)[8]="percent.change.DOC"
names(total.df)[9]="abs.change.TN"
names(total.df)[10]="percent.change.TN"
names(total.df)[11]="SUVA"
names(total.df)[12]="Spectral.slope"
names(total.df)[13]="HIX"
names(total.df)[14]="Fr"
names(total.df)[15]="FI"
total.df$mass.remaining.DOC=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.DOC)
total.df$mass.remaining.TN=100-abs(total.df[total.df$effect=="UV and microbes",]$percent.change.TN)
total_t1.df=subset(total.df, Time.point %in% c("T1"))
total.df=total.df%>%
mutate_if(is.character,as.factor)
total_t1.df=total_t1.df%>%
mutate_if(is.character,as.factor)
#merge degra.df time points into one dataframe
degra.df=rbind(photo.df, microbe.df, total.df)
degra.df=pivot_longer(degra.df, cols=c("effect"), names_to="effect", values_to="Type")
degra_t1.df=rbind(photo_t1.df, microbe_t1.df, total_t1.df)
degra_t1.df=pivot_longer(degra_t1.df, cols=c("effect"), names_to="effect", values_to="Type")
#make sure we have factors
degra.df=degra.df%>%
mutate_if(is.character,as.factor)
degra_t1.df=degra_t1.df%>%
mutate_if(is.character,as.factor)
#seperate degra.df for burned and unburned tanks
degra_t1.burned.df=degra_t1.df[degra_t1.df$Treatment=="burned",]
degra_t1.unburned.df=degra_t1.df[degra_t1.df$Treatment=="unburned",]UVBIO analysis
We can model how UVBIO changes across the biomass gradient on Day 10.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for UVBIO.
library(mgcv)
###Just UV
m1.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m2.photo.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m3.photo.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m1.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m2.photo.abs=gam(abs.change.DOC ~ Treatment +s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
m3.photo.abs=gam(abs.change.DOC ~ s(plant.mass..g), data=photo.df, method="REML", select= T, subset= Time.point=="T1")
photo.AIC=AIC(m1.photo.perc, m2.photo.perc, m3.photo.perc)
photo.AIC=AIC(m1.photo.abs, m2.photo.abs, m3.photo.abs)
photo.AIC
summary(m1.photo.perc)
photo.anova=anova(m1.photo.perc)
photo.anova
gam.check(m3.photo.perc, rep=1000)
draw(m3.photo.perc)
concrvity(m3.photo.perc)
par(mfrow = c(2, 2))
plot(m3.photo.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
photo.mod.plot<-
plot_smooths(
model = m1.photo.perc,
series = plant.mass..g,
comparison = Treatment
) +
geom_point(data=photo.df[(photo.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
geom_hline(yintercept=0, linetype="longdash", color = "gray")+
ggtitle("UV-only") +
ylab("Change DOC (%)") +
xlab("plant material (g)") +
coord_cartesian(ylim = c(-30,10))+
Fig.formatting
photo.mod.plot
###Just microbes
m1.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
m2.microbe.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
m3.microbe.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=microbe.df, method="REML", select= T, subset= Time.point=="T1")
microbe.AIC=AIC(m1.microbe.perc, m2.microbe.perc, m3.microbe.perc)
microbe.AIC
summary(m3.microbe.perc)
microbe.anova=anova(m1.microbe.perc)
gam.check(m3.microbe.perc, rep=1000)
draw(m3.microbe.perc)
concrvity(m3.microbe.perc)
par(mfrow = c(2, 2))
plot(m3.microbe.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
microbe.mod.plot<-
plot_smooths(
model = m3.microbe.perc,
series = plant.mass..g
) +
geom_point(data=microbe.df[(microbe.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("microbes-only") +
geom_hline(yintercept=0, linetype="longdash", color = "gray") +
ylab("Change DOC (%)") +
coord_cartesian(ylim = c(-30,10))+
xlab("plant material (g)") +
Fig.formatting
coord_cartesian(ylim = c(-30,0))
###Just total decomp
m1.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g, by= Treatment), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m2.total.perc=gam(percent.change.DOC ~ Treatment +s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m3.total.perc=gam(percent.change.DOC ~ s(plant.mass..g), data=total.df, method="REML", select= T, subset= Time.point=="T1")
m4.total.perc=gam(percent.change.DOC ~ s(plant.mass..g, m=c(2,0)), data=total.df, method="REML", select= T, subset= Time.point=="T1")
total.AIC=AIC(m1.total.perc, m2.total.perc, m3.total.perc)
total.AIC
summary(m3.total.perc)
total.anova=anova(m1.total.perc)
gam.check(m3.total.perc, rep=1000)
draw(m3.total.perc)
concrvity(m3.total.perc)
par(mfrow = c(2, 2))
plot(m3.total.perc, all.terms = TRUE, page=1)
#plot for the model output on rawdata
total.mod.plot<-
plot_smooths(
model = m3.total.perc,
series = plant.mass..g
) +
geom_point(data=total.df[(total.df$Time.point=="T1"),],
aes(x=plant.mass..g, y=percent.change.DOC, color=Treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("UV and microbes") +
geom_hline(yintercept=0, linetype="longdash", color = "gray") +
ylab("Change DOC (%)") +
coord_cartesian(ylim = c(-30,10))+
xlab("plant material (g)") +
Fig.formatting
##combine AIC
mod.degra.rep<-rep(c("UV-only: Treatment with by-factor smooth",
"UV-only: Treatment with global smooth",
"UV-only: global smooth",
"Microbes-only: Treatment with by-factor smooth",
"Microbes-only: Treatment with global smooth",
"Microbes-only: global smooth",
"UV and microbes: Treatment with by-factor smooth",
"UV and microbes: Treatment with global smooth",
"UV and microbes: global smooth"))
mod.degra.column<- data.frame(mod.degra.rep)
AIC.degra<-bind_rows(photo.AIC, microbe.AIC, total.AIC)
AIC.degra.mod<-as.data.frame(cbind(mod.degra.column, AIC.degra))
AIC.degra.mod
names(AIC.degra.mod)[1]="Model"
AIC.degra.mod <- AIC.degra.mod %>%
# Creating an empty column:
add_column(Variable= c("-"), .before="Model")
write.csv(AIC.degra.mod, "output/AIC models/AIC.degra.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.AIC.pdf", width = 6.5, height = 2.8) # Export PDF
grid.table(AIC.degra.mod, rows=NULL)
dev.off()
library(effectsize)
r_squared_microbes <- summary(m1.microbe.perc)$r.sq
r_squared_total <- summary(m1.total.perc)$r.sq
r_squared_microbes
r_squared_total
UVBIO output
This generates a table with GAM summaries.
####Parametric output
degra.para.sums=as.data.frame(rbind(photo.anova$pTerms.table, microbe.anova$pTerms.table, total.anova$pTerms.table))
names(degra.para.sums)[1]="df/edf"
#create empty column for ref.df
degra.para.sums <- degra.para.sums %>%
# Creating an empty column:
add_column(Model= c("UV-only Model: Treatment","Microbes-only Model: Treatment","UV and microbes Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
degra.smooth.sums=as.data.frame(rbind(photo.anova$s.table, microbe.anova$s.table, total.anova$s.table))
names(degra.smooth.sums)[1]="df/edf"
#create empty column for ref.df
degra.smooth.sums <- degra.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("UV-only Model: burned smooth","UV-only Model: unburned smooth","Microbes-only Model: burned smooth","Microbes-only Model: unburned smooth","UV and microbes Model: burned smooth","UV and microbes Model: unburned smooth"), .before="df/edf")
#merge the dataframes
degra.mod.allsums=rbind(degra.para.sums,degra.smooth.sums)
degra.mod.allsums <- degra.mod.allsums %>%
# Creating an empty column:
add_column(Variable= c("-"), .before="Model")
write.csv(degra.mod.allsums,"/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/degra.mod.allsums.pdf", width = 8.5, height = 2.8) # Export PDF
grid.table(degra.mod.allsums, rows=NULL)
dev.off()UVBIO model plot
degra.mod.plots = plot_grid(
photo.mod.plot + theme(legend.position = "none"),
microbe.mod.plot + theme(legend.position = "none"),
total.mod.plot + theme(legend.position = "none" ),
extract.legend,
rel_widths = c(8,8,8,3), nrow = 1)
degra.mod.plots
Figure. Dissolved organic carbon (DOC) loss as percent change for burned (red) and unburned (green) treatments across plant biomass (g). Each point represents a different pond. Plots with a single smooth (global) indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/degra.mod.plots.pdf", height=4, width=10)Import sage and willow decomposition data
Plant material for the starting material (sage or willow). This is useful in determining how fire impacted decomposition of particulate detritus.
- First, run GAMs to see what is happening and where differences
lie.
- We will then make a plot for sage and willow within and between treatments.
#load file
plant.df=read.csv("/Users/Cody1/Desktop/GitHub/Pyromania2/data/decomp.weights.csv")
plant.df=plant.df[-c(61:82), ]
plant.df=plant.df[-c(6, 8,9)]
#make sure we have factors
plant.df=plant.df%>%
mutate_if(is.character,as.factor)
#add propotion column so it can be arcsine transformed
plant.df$perc.loss=(100*(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g-plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$dry.mass..g.bag.corrected)/(plant.df[plant.df$plant.type=="sage" | plant.df$plant.type=="willow",]$plant.specific.mass..g))
#remove NaN
plant.df[is.na(plant.df)] <- 0
#calculate mass remaining
plant.df$mass.remaining=100-plant.df$perc.loss
hist(plant.df$mass.remaining)
#add delta
plant.df$delta.loss= plant.df$perc.loss/100
#arcsine
plant.df$delta.trans=asin(sqrt(plant.df$delta.loss))
hist(plant.df$delta.trans)
#exclude control tank outlier
plant.df2=data.frame(plant.df[plant.df$plant.specific.mass..g > 0,])
###Sage
m1.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m2.sage.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m3.sage.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m4.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
hist(m3.sage.gam$residuals)
m5.sage.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
m6.sage.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="sage")
sage.AIC=AIC(m4.sage.gam, m5.sage.gam, m6.sage.gam)
sage.AIC
sage.anova=anova(m4.sage.gam)
sage.anova
gam.check(m6.sage.gam, rep=1000)
mean(plant.df[plant.df$plant.type=="sage",]$delta.trans)
# global smooth
# model predictions
sage.diff.mod<-plot_difference(
m1.sage.gam,
series = plant.specific.mass..g,
difference = list(treatment = c("burned", "unburned"))
)
sage.diff.mod
summary(m4.sage.gam)
anova.gam(m3.sage.gam)
gam.check(m4.sage.gam, rep=1000)
draw(m3.sage.gam)
concrvity(m3.sage.gam)
par(mfrow = c(2, 2))
plot(m3.sage.gam, all.terms = TRUE, page=1)
qq_plot(m3.sage.gam)
#plot for the model output on rawdata
sage.mod.plot<-
plot_smooths(
model = m4.sage.gam,
series = plant.specific.mass..g,
comparison= treatment
) +
geom_point(data=plant.df2[(plant.df2$plant.type=="sage"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Sage") +
ylab("Mass loss (%)") +
coord_cartesian(ylim = c(0,80))+
xlab("Plant material (g)") +
Fig.formatting
sage.mod.plot
### Willow
m1.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m2.willow.gam=gam(delta.trans ~ treatment +s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m3.willow.gam=gam(delta.trans ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m4.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g, by= treatment), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
m5.willow.gam=gam(perc.loss ~ treatment +s(plant.specific.mass..g), data=plant.df, method="REML", select= T, subset= plant.type=="willow")
m6.willow.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, method="REML", select= T, subset= plant.type=="willow")
willow.AIC=AIC( m4.willow.gam,m5.willow.gam,m6.willow.gam)
willow.AIC
# treatment with global smooth
summary(m4.willow.gam)
willow.anova=anova(m4.willow.gam)
willow.anova
gam.check(m3.willow.gam, rep=1000)
draw(m2.willow.gam)
concrvity(m2.willow.gam)
par(mfrow = c(2, 2))
plot(m2.willow.gam, all.terms = TRUE, page=1)
# model predictions
willow.diff.mod<-plot_difference(
m1.willow.gam,
series = plant.specific.mass..g,
difference = list(treatment = c("burned", "unburned"))
)
willow.diff.mod
#plot for the model output on rawdata
willow.mod.plot<-
plot_smooths(
model = m4.willow.gam,
series = plant.specific.mass..g,
comparison = treatment
) +
geom_point(data=plant.df2[(plant.df2$plant.type=="willow"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=treatment)) +
scale_color_manual(values = c("brown1", "mediumseagreen")) +
scale_fill_manual(values = c("brown1", "mediumseagreen")) +
ggtitle("Willow") +
ylab("Mass loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
willow.mod.plot
##combine AIC
mod.decomp.rep<-rep(c("Sage: Treatment with by-factor smooth",
"Sage: Treatment with global smooth",
"Sage: global smooth",
"Willow: Treatment with by-factor smooth",
"Willow: Treatment with global smooth",
"Willow: global smooth"))
mod.decomp.column<- data.frame(mod.decomp.rep)
AIC.decomp<-bind_rows(sage.AIC, willow.AIC)
AIC.decomp.mod<-as.data.frame(cbind(mod.decomp.column, AIC.decomp))
names(AIC.decomp.mod)[1]="Model"
write.csv(AIC.decomp.mod, "output/AIC models/AIC.decomp.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.AIC.pdf", width = 5, height = 2) # Export PDF
grid.table(AIC.decomp.mod, rows=NULL)
dev.off()
Sage and willow output
This generates a table with GAM summaries.
####Parametric output
decomp.para.sums=as.data.frame(rbind(sage.anova$pTerms.table, willow.anova$pTerms.table))
names(decomp.para.sums)[1]="df/edf"
#create empty column for ref.df
decomp.para.sums <- decomp.para.sums %>%
# Creating an empty column:
add_column(Model= c("Sage Model: Treatment","Willow Model: Treatment"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
decomp.smooth.sums=as.data.frame(rbind(sage.anova$s.table, willow.anova$s.table))
names(decomp.smooth.sums)[1]="df/edf"
#create empty column for ref.df
decomp.smooth.sums <- decomp.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Sage Model: burned smooth","Sage Model: unburned smooth","Willow Model: burned smooth","Willow Model: unburned smooth"), .before="df/edf")
#merge the dataframes
decomp.mod.allsums=rbind(decomp.para.sums,decomp.smooth.sums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.mod.allsums.pdf", width = 7, height = 2) # Export PDF
grid.table(decomp.mod.allsums, rows=NULL)
dev.off()
Sage and willow within treatment analysis
Above, we examined differences between plant types. We can take a closer look to see how sage and willow decomposition differed within fire treatment (i.e., differences between burned sage and willow, and unburned sage and willow).
#compare mass loss between sage and willow by fire treatment---burned
m1.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
m2.compare.burned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
m3.compare.burned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "burned")
compare.burned.AIC=AIC(m1.compare.burned.gam, m2.compare.burned.gam, m3.compare.burned.gam)
compare.burned.anova=anova(m1.compare.burned.gam)
compare.burned.anova
gam.check(m1.compare.burned.gam, rep=1000)
# model predictions
compare.burned.diff.mod<-plot_difference(
m1.compare.burned.gam,
series = plant.specific.mass..g,
difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod
#plot for the model output on rawdata
compare.burned.mod.plot<-
plot_smooths(
model = m2.compare.burned.gam,
series = plant.specific.mass..g,
comparison = plant.type
) +
geom_point(data=plant.df2[(plant.df2$treatment=="burned"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
ggtitle("Burned") +
geom_line(aes(linetype= plant.type))+
ylab("Mass loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
compare.burned.mod.plot
#compare mass loss between sage and willow by fire treatment---unburned
m1.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g, by=plant.type), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
m2.compare.unburned.gam=gam(perc.loss ~ plant.type + s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
m3.compare.unburned.gam=gam(perc.loss ~ s(plant.specific.mass..g), data=plant.df2, methods= "REML", select= T, subset= treatment == "unburned")
compare.unburned.AIC=AIC(m1.compare.unburned.gam, m2.compare.unburned.gam, m3.compare.unburned.gam)
compare.unburned.anova=anova(m1.compare.unburned.gam)
compare.unburned.anova
summary(m3.compare.unburned.gam)
gam.check(m3.compare.unburned.gam, rep=1000)
mean(plant.df2[plant.df2$plant.type=="sage",]$perc.loss)
mean(plant.df2[plant.df2$plant.type=="willow",]$perc.loss)
# model predictions
compare.burned.diff.mod<-plot_difference(
m1.compare.burned.gam,
series = plant.specific.mass..g,
difference = list(plant.type = c("sage", "willow"))
)
compare.burned.diff.mod
#plot for the model output on rawdata
compare.unburned.mod.plot<-
plot_smooths(
model = m1.compare.unburned.gam,
series = plant.specific.mass..g,
comparison = plant.type
) +
geom_point(data=plant.df2[(plant.df2$treatment=="unburned"),],
aes(x=plant.specific.mass..g, y=perc.loss, color=plant.type)) +
scale_color_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
scale_fill_manual(values = c("mediumpurple4", "darkgoldenrod4")) +
ggtitle("Unburned") +
ylab("Mass Loss (%)") +
xlab("Plant material (g)") +
coord_cartesian(ylim = c(0,80))+
Fig.formatting
compare.unburned.mod.plot
##combine AIC
mod.comapre.decomp.rep<-rep(c("Burned: plant type with by-factor smooth",
"Burned: plant type with global smooth",
"Burned: global smooth",
"Unburned: plant type with by-factor smooth",
"Unburned: plant type with global smooth",
"Unburned: global smooth"))
mod.compare.decomp.column<- data.frame(mod.comapre.decomp.rep)
AIC.compare.decomp<-bind_rows(compare.burned.AIC, compare.unburned.AIC)
AIC.compare.decomp.mod<-as.data.frame(cbind(mod.compare.decomp.column, AIC.compare.decomp))
names(AIC.compare.decomp.mod)[1]="Model"
write.csv(AIC.compare.decomp.mod, "output/AIC models/AIC.comapre.decomp.csv")
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.AIC.pdf", width = 5.3, height = 2) # Export PDF
grid.table(AIC.compare.decomp.mod, rows=NULL)
dev.off()
#### generate plots
decomp.compare.mod.plot<-plot_grid(
compare.burned.mod.plot+ theme(legend.position = "none"),
compare.unburned.mod.plot+ theme(legend.position = "none"),
ncol=2)
ggsave("figures/decomp.comapre.mod.plot.pdf", height=4, width=12)
decomp.mod.plot<-plot_grid(
sage.mod.plot+ theme(legend.position = "none"),
willow.mod.plot+ theme(legend.position = "none"),
rel_widths = c(8,8), ncol=2)
ggsave("figures/decomp.mod.plot.pdf", height=4, width=12)
Sage and willow within treatment output
This generates a table with GAM summaries.
####Parametric output
decomp.compare.para.sums=as.data.frame(rbind(compare.burned.anova$pTerms.table, compare.unburned.anova$pTerms.table))
names(decomp.compare.para.sums)[1]="df/edf"
#create empty column for ref.df
decomp.compare.para.sums <- decomp.compare.para.sums %>%
# Creating an empty column:
add_column(Model= c("Burned: Plant type","Unburned: Plant type"), .before="df/edf")%>%
add_column(Ref.df = "-", .after="df/edf")
#####Smooth output
decomp.compare.smooth.sums=as.data.frame(rbind(compare.burned.anova$s.table, compare.unburned.anova$s.table))
names(decomp.compare.smooth.sums)[1]="df/edf"
#create empty column for ref.df
decomp.compare.smooth.sums <- decomp.compare.smooth.sums %>%
# Creating an empty column:
add_column(Model= c("Burned: sage smooth","Burned: willow smooth","Unburned Model: sage smooth","Unburned Model: willow smooth"), .before="df/edf")
#merge the dataframes
decomp.compare.mod.allsums=rbind(decomp.compare.para.sums,decomp.compare.smooth.sums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/decomp.compare.mod.allsums.pdf", width = 6.2, height = 2) # Export PDF
grid.table(decomp.compare.mod.allsums, rows=NULL)
dev.off()
Sage and willow combined output
This combines to outputs into a single table with GAM summaries.
library(gridExtra)
#merge all the dataframes
all.decomp.mod.sums=rbind(decomp.mod.allsums,decomp.compare.mod.allsums)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.sums.pdf", width = 7, height = 3.7) # Export PDF
grid.table(all.decomp.mod.sums, rows=NULL)
dev.off()
#merge all the dataframes
all.decomp.mod.AIC=rbind(AIC.decomp.mod, AIC.compare.decomp.mod)
#convert dataframe into a table for export
pdf("/Users/Cody1/Desktop/GitHub/Pyromania2/output/Codys.tables/all.decomp.mod.AIC.pdf", width = 5.3, height = 3.7) # Export PDF
grid.table(all.decomp.mod.AIC, rows=NULL)
dev.off()Sage and willow combined plot
We can combine the two decomposition plots into one so that we now have decomposition between and within fire treatment.
##### combine all plots
all.decomp.mod.plot<-plot_grid(
decomp.mod.plot+ theme(legend.position = "none"),
decomp.compare.mod.plot+ theme(legend.position = "none"),
nrow=2)
all.decomp.mod.plot
Figure. Sage and willow percent mass loss across the plant biomass gradient at the end of 89 days. Mass loss for sage and willow between burned (red) and unburned (green) treatments. Mass loss between sage (purple) and willow (brown) within treatments. Each point represents a different pond. Plots with a single global smooth indicate that a model with one fit line for the burned and unburned treatments best describes the data, while separate smooths indicate that a model with different fits for the two treatments showed best agreement.
ggsave("figures/all.decomp.mod.plot.pdf", height=8, width=8)